Extended Rate Constants Distribution (RCD) Model for Sorption in Heterogeneous Systems: 4. Kinetics of Metal Ions Sorption in the Presence of Complexing Agents—Application to Cu(II) Sorption on Polyethyleneimine Cryogel from Acetate and Tartrate Solutions

Here, we report a new version of the extended Rate Constants Distribution (RCD) model for metal ion sorption, which includes complex-formation equilibria. With the RCD-complex model, one can predict sorbent performance in the presence of complexing agents using data on metal ion sorption from ligand-free solutions and a set of coefficients for sorption rate constants of different ionic species. The RCD-complex model was applied to breakthrough curves of Cu(II) sorption from acetate and tartrate solutions on polyethyleneimine (PEI) monolith cryogel at different flow rates and ionic speciation. We have shown that, despite the lower stability of Cu(II)-acetate complex, at high flow rates, acetate has a more pronounced negative effect on sorption kinetics than tartrate. The RCD model was successfully used to predict the shape of the breakthrough curves at an arbitrary acetate concentration but failed to predict Cu(II) sorption from tartrate solutions in a broad range of ligand concentrations. Since a twofold increase in sorption capacity was observed at low tartrate concentrations, the latter fact was related to an alteration in the sorption mechanism of Cu(II)-ions, which depended on Cu(II) ionic speciation. The obtained results emphasize the importance of information about sorption kinetics of different ionic forms for the optimization of sorption filter performance in the presence of complexing agents.


Introduction
The problem of heavy metal removal from natural and waste water remains in focus, despite technological achievements and the availability of a broad selection of different types of sorbents [1][2][3][4][5][6]. Aside from affinity, selectivity, and high sorption capacity, the kinetic characteristics of sorbents are of utmost importance as they determine the productivity of both industrial water treatment facilities and small-size point-of-use (POU) filters. Demand for materials, which can efficiently remove metal ions at high flow rates in fixed-bed applications has resulted in a sharp growth in demand for highly porous sorbents that are based on water-soluble synthetic and natural polymers [4][5][6][7][8][9][10][11]. However, when diffusion limitations are eliminated due to the well-developed porous structure, chemisorption (surface chemical reaction) can become a rate-limiting stage of adsorption. In this case, even in single metal solutions, metal speciation can significantly affect not only equilibrium sorption parameters but also sorption kinetics via an alteration of the sorption mechanism and a difference in the sorption/desorption rate constants for different species.
Ligands can alter both metal speciation and solid surface properties, leading to negative or positive cooperative effects on adsorption, depending on adsorption mechanisms, heterogeneous sorbents via simultaneous processing of several experimental curves of sorption kinetics in batch or in fixed-bed applications. This single RCD function, which describes the fill set of experimental curves, contains all information about affinity, quantity, and distribution of the sorption sites in the space of constants of sorption and desorption rates for heterogeneous sorbent and can be used to predict distribution of the adsorbate on sorption centers for different starting conditions (solid:liquid ratio or column geometry, adsorbate concentration, and flow rate) at any time of the sorption process.
Recently, we reported the fabrication of polyethylenimine (PEI) monolith supermacroporous cryogels with high efficiency of heavy metal ions under dynamic conditions at flow rates up to 160 bed volume (BV)/h [35]. Application of the RCD model to Cu(II), Ni(II), Zn(II), and Cd(II) sorption on PEI cryogel in batch and fixed-bed experiments [27,29] allowed identification of "fast" and "slow" sorption centers. This model also provided an explanation for the preferential adsorption of one or another ion from the mixture in a fixedbed application, depending on the experimental conditions such as flow rate and metal ion concentrations [29]. Moreover, RCD functions for Cu(II), Zn(II), and Cd(II) sorption on PEI cryogels obtained from the batch data were successfully used for the first time to predict the rate dependence of a breakthrough point position and shape of the breakthrough curves in a fixed-bed application for a broad range of adsorbate concentrations and flow rates [29].
However, in the presence of complexing agents, metal ion speciation in solutions differs depending on metal/ligand ratio and pH; thus, RCD functions obtained for metal cations sorption in ligand-free solutions cannot be used to model sorption kinetics in the presence of ligands. Obviously, another RCD function can be obtained for the sorption of a metal complex, but its application will be limited to the systems containing only one ionic form of the metal and will not be applicable to the broad range of metal/ligand ratios and mixture of ionic species. This problem can be resolved if we extend the RCD model with complex-formation equilibria and find constants of sorption/desorption rates for all possible ionic forms. Mathematically, it can be more easily performed for the systems for which binding constants of a metal ion to the sorbent functional groups will significantly exceed the binding constant to a ligand. In other words, the structure of the metal-sorbent complex will be the same after the sorption from ligand-free solution and solutions containing different metal-ligand species. So, if one can somehow modify RCD function for the metal ion sorption from ligand-free solution to take into account complexformation equilibria, one will be able to perform predictive modeling of the breakthrough curves of metal sorption at various flow rates and metal/ligand ratios.
Such systems, containing ligands, which form relatively weak complexes with transition metal ions, are common in industry and in the environment. One example is lowmolecular-weight carboxylic acids, which are extensively used for soil remediation, leaching of metal ores [36], electroplating [37], and circuit-board printing, and thus shall be considered as co-contaminants in metal ion removal by sorption methods. The effect of their presence on the kinetics of metal ion removal in a fixed-bed application at high flow rates is important for wastewater treatment and an estimation of the sorption filters' productivity.
Thus, in this work, which continues a series of papers on RCD model development and applications [27][28][29], we investigated how the presence of acetates and tartrates affects Cu(II) sorption in a fixed-bed application on a monolith supermacroporous PEI cryogel at different flow rates and metal/ligand ratios. Moreover, we verified whether the possibility of the RCD model extension to predict kinetics of the metal ion sorption in the presence of complexing agents is dependent on the sorption mechanism.

Brief Introduction to RCD Model
The Extended Rate Constants Distribution (RCD) model for sorption in heterogeneous systems belongs to the pool of models, which assumes the presence of a continuum of different sorption sites [38][39][40][41]. The background and advantages of this model were verified and discussed in our earlier works [27][28][29]. The following assumptions were made in the RCD model: (1) the flow of adsorbate from the bulk to the sorbent is proportional to the  adsorbate concentration in the solution and to the surface area with vacant sorption sites; (2) the flow of adsorbate from the surface to the bulk is proportional to the surface area with occupied sorption sites; and (3) the specific surface area occupied with sorption sites of the certain type is proportional to the content of such sites in the sorbent.
The differential mass balance equation for a fixed-bed column is given in [42]: where u is the superficial velocity, C is the adsorbate concentration in the solution, z is the axial coordinate, τ is the time, ε is the column void fraction, ρ is the adsorbent density, Q is the adsorbate content in the sorbent, and D L is the axial dispersion coefficient, which was set here to zero as in most other works on sorption dynamics [42,43]. Any suitable rate expression for ∂Q/∂τ can be used to complete the fixed-bed model. Because, in most cases, metal ion adsorption is well described by Langmuir equation, in the RCD model we have used Langmuir adsorption kinetic Equation (2): where k s is the second-order sorption rate constant, k d is the first-order desorption rate constant, C-the adsorbate concentration, Q max is the total sorption capacity, and Q(τ) is the total content of the adsorbate in the sorbent at time τ.
To describe kinetics of sorption in the RCD model, we have introduced the density function q(k s , k d , τ), which shows the distribution of the adsorbate on sorption sites in the space of the rate constants (RC) of sorption (k s ) and desorption (k d ) at any time point (τ), and the density function q max (k s , k d ), which shows the maximal content of the adsorbate (sorption capacity) for the certain type of sorption sites (k s , k d ) at full saturation. Q max is the total sorption capacity. Using these functions, Equation (2) will transform to Equation (3): After addition of the material balance, Equation (4): where Q 0 is adsorbate content in the sorbent, C 0 is the adsorbate concentration in the solution in the initial time point, and V sp is the specific solution volume. For the heterogeneous sorbents, we can write the following integral equations, where Equation (6) describes the Langmuir-type sorption isotherm on a heterogeneous sorbent: where C e and Q e are the equilibrium concentration of the adsorbate in the solution and the content of the adsorbate in the sorbent, respectively. Transformation of the systems of integro-differential Equations (3)- (6) to the form suitable for numerical calculations and processing the experimental data are described in [27].

RCD Model for Sorption in the Presence of Complexing Agents (RCD-Complex Model)
Let us consider metal ion (Me) sorption on heterogeneous sorbent in the presence of complex-forming ligands (L). Assuming that n-maximal coordination number of metal ions, the following metal ionic forms will be present in solution in addition to the free metal cations: MeL, MeL 2 , . . . , MeL n . The following equilibria can be written without taking into account ions charge (this is valid for the case when L-is not hydroxyl ion): k 1 . . . k n -the stepwise complex instability constants. Material balance equation for metal is: where v-the solution volume per mass unit of the sorbent; q i -the content of metal on the i-the sorption center; m-the number of sorption centers; and N 0 Me -the initial number of moles of metal in solution in sorbent. The equation of material balance with respect to the ligand is: N 0 L -the initial number of moles of a ligand; C L -the current concentration of free ligands.
If we assume momentary establishment of equilibrium concentrations of ionic forms in the solution, then to calculate equilibrium concentrations of ionic forms in the initial solution (without the sorbent), the equations of the material balance must be supplemented with the law of acting masses for all ionic forms: Thus, we obtain the following system of nonlinear equations: Equation (11) was solved by the Newton-Raphson method to ensure the positive sign of concentrations of all ionic forms. Through solving Equation (11), one obtains the distribution in the forms of existence of metal ions for the initial time point. At the other time points, Equation (11) in differential forms was included in the system of differential equations of sorption kinetics and considered during integration: Let us assume that (1) the final form of binding of any complex metal forms with the sorbent surface-Me; (2) all the complex ionic forms Me are sorbed (with different desorption rates until 0); and (3) the constants of sorption rates for metal complex forms are related through functional dependencies with sorption/desorption rates of the ion Me (if we know the function of distribution of metal ions without the complex ion). Therefore, the sorption equation (just like Equation (3) in Section 2.1 of the Langmuir sorption kinetics in the absence of the complexing agent) can be written as: where K s,MeL i (pk s,Me , pk d,Me )-the function of dependence of the rate of sorption of the ionic form MeL i on logarithms of constants of sorption/desorption rates of the ion Me; q Me (k s,Me , k d,Me )-the function of distribution of the content of the ion Me over sorption centers; and C Me (τ), C MeL i (τ)-the concentration of ion and ionic forms of Me in solution in the time τ. The appearance of the functions K s,MeL i (pk s,Me , pk d,Me ) is an unknown a priori, but they can be expanded into a Taylor series around the point (pk s,Me , pk d,Me ): K s (pk s , pk d ) = a 00 + a 10 pk s + a 01 pk d + a 20 pk 2 s + a 11 pk s pk d + a 02 pk 2 d + . . .
where a ij -the coefficient of a two-dimensional Taylor series by the term k i s k j d (see Appendix A for details).
To sum up, if one knows the function of distribution of sorption centers for the ion Me, then to describe sorption of metal complex forms, one must determine the coefficients a ij for every complex form. The latter can be easily realized through minimization of the functional below (taking into account the solution regularization): where a MeL i -the vector of unknown coefficients in the expansion (14) for the ionic form MeL i ; M, N i -the number of experimental kinetic curves for sorption of complex forms of Me and the number of experimental points on the curve i; C ij,calc (a MeL ,a MeL 2 , . . .)-the calculated value of the metal ion concentration in the point i on the curve j for specific values of the vector components a ij ; Ω(a MeL ,a MeL 2 , . . .)-the solution stabilizer; λ-the regularization parameter; a MeL k − a 0,MeL k ¯the Euclid form of the difference between the "zero" and current vector a MeL k (i.e., it is preferable that the optimal vector of expansion coefficients would not differ significantly from the "zero" one); and n L -the maximal degree of complex binding of the metal ion.
To sum up, in order to describe the metal sorption in the presence of the complexing agent, it is first necessary (i) to determine the function of density of sorption centers for sorption of pure Me; (ii) to obtain the number of sorption kinetic curves from metal-chelate solutions, preferably, differing in metal-ligand concentrations; (iii) to find an adequate set of coefficients in expansion (14) via the functional minimization (15). The initial values of the vector a ij components are determined by the conventional method (the constants of rates for all the complex forms are equal to those of metal ions): In other words, the initial values of constants of rated for all the ionic forms are equal to those of metal ions.

Cu(II) Soprtion on PEI Cryogel in the Presence of Acetate and Tartrate
Fabrication and characterization of PEI cryogels using diglycidyl ether of 1,4-butanediol (DGEBD) as a cross-linker was reported earlier in [35]. At a molar ratio DGEBD: PEI of 1:4, we obtained a highly permeable monolith cryogel with a swelling degree of 2200% and a pore size of 128 ± 30 µm, which supported a liquid flow rate up to 450 BV/h [35]. Efficient mass transfer through the channels of interconnected macropores under dynamic conditions and very low thickness (6.1 ± 2.6 µm) of the polymeric walls assured very high rates of metal ion sorption in this cryogel via coordination mechanism [29]. However, in case of electrostatic interactions between positively charged PEI and anionic adsorbates (alizarin red dye and humic acid), an increase in the surface coverage reduced attractive forces between the adsorbate, so the sorption rate slowed down and breakthrough curves featured an asymmetric shape with a bend after the breakpoint and significant difference between effective dynamic sorption capacities at moderate flow rates.
To investigate how the presence of acetic and tartaric acids as complexing agents affect ability of PEI cryogel to remove Cu(II) ions under dynamic conditions, we have designed several model solutions with different compositions of metal ionic forms using Phreeqc Interactive 3.7.3-15968 software (Tables 1 and 2). Taking into account that wastewaters can contain a much higher excess of the complexing agents than were earlier used for studying Cu(II) sorption on PEI and polyamine resin in the presence of citric acid and EDTA (metal:ligand 1:2) [16,18,19], acetate and tartrate concentrations were varied over a broad range.  First, breakthrough curves of Cu(II) sorption from 0.1 M acetate and tartrate solutions were obtained at several flow rates from 8 to 163 BV/h and compared with those for Cu(II) sorption from water ( Figure 1). Since the intraparticle diffusion in monolith PEI cryogel is not the rate-limiting stage of sorption [29,35], the difference in the breakthrough curve shape (slop and breakthrough point position) in water and ligand-containing solutions will reflect the difference in sorption rate constants of Cu(II) ionic species and a possible alteration of sorption mechanism. At flow rates of 41 BV/h and 163 BV/h, the effective sorption capacities of PEI for Cu(II) from 0.1 M acetate were only 80% and 50% of those in water, respectively. This shows that Cu(II) speciation in 0.1 M acetate solution significantly affects sorption kinetics despite the low thermodynamic stability of Cu(II)-acetate complex. From a practical point, this means that the presence of acetate will decrease filter productivity if the target of sorption process is to meet requirements for Cu content in discharge waters.  (Table  2), which can be adsorbed by PEI via strong electrostatic interactions without a ligand- Both acetate and tartrate Cu(II) complexes are significantly weaker in comparison with Cu(II)-citrate and Cu(II)-EDTA complexes (cumulative stability constant 14.2 and 18.7, respectively [44]). So, one can expect that, at equilibrium, Cu(II) ions will be bound by PEI cryogel with the release of the ligands to the solution. Considering the higher stability of Cu(II) tartrate complexes, tartrate could have a more profound effect on Cu(II) sorption than acetate. However, the effective dynamic sorption capacities for Cu(II) in 0.1 M tartrate were similar or even higher than those obtained in water, while the efficiency of Cu(II) uptake from 0.1 M acetate solution significantly dropped with an increasing flow rate ( Figure 1, Table 3). The effective dynamic capacities were calculated for the breakthrough point of 2 mg/L that corresponds to the World Health Organization guideline value for copper in drinking water [45]. tartrate (3.87 over 1.73 mmol/g), even at a higher flow rate of 130 BV/h (Table 3). These data are not in agreement with conclusion made for Cu(II) sorption from citrate solutions on polyaminated resin about the higher sorption rate constant for the neutral CuHL 0 complex over anionic forms [16]. Regretfully, experimental data in [16] were obtained in batch and were not verified in fixed-bed models at different flow rates and metal/ligand ratios that could prove one or another hypothesis. At the same time, and similarly to the works [16,18,19,46] on Cu(II) sorption on polyamine resins in the presence of citrate, we observed a drastic increase in maximal sorption capacity for Cu(II) in the presence of a low excess of tartrate ( Table 3). The maximal value of 5.61 mmol/g is in very good agreement with the sorption capacity of 5.5 mmol/g toward anionic dye earlier determined for this cryogel [35]. It should be mentioned that homogeneous distribution of Cu(II) in monolith cryogel after the sorption from both acetate and tartrate solutions was confirmed by SEM ( Figure S1, Supplementary information). Thus, the difference in sorption capacities and rate-dependence of the shape of the breakthrough curves can be attributed to the different mechanisms of Cu(II) sorption in the presence of acetate and tartrate. At flow rates of 41 BV/h and 163 BV/h, the effective sorption capacities of PEI for Cu(II) from 0.1 M acetate were only 80% and 50% of those in water, respectively. This shows that Cu(II) speciation in 0.1 M acetate solution significantly affects sorption kinetics despite the low thermodynamic stability of Cu(II)-acetate complex. From a practical point, this means that the presence of acetate will decrease filter productivity if the target of sorption process is to meet requirements for Cu content in discharge waters.
In 0.1 M tartrate solution, Cu(II) exists predominantly in anionic form [CuL 4 ] 6− ( Table 2), which can be adsorbed by PEI via strong electrostatic interactions without a ligand-exchange reaction and the release of tartrate to the solution. This is a fast process and at the highest flow rate of 163 BV/h, effective dynamic sorption capacity for Cu(II) in tartrate is only 22% lower than that in water (Table 3). When tartrate concentration was decreased to a Cu(II):L ratio of 1:2, Cu(II) ionic speciation shifted toward a higher contribution (41.9%) of neutral form (0.00312 M tartrate, Table 2), and the increase in flow rate had a more profound effect of effective dynamic sorption capacity, which was threefold lower at 84 BV/h than at 8 BV/h ( Figure 2a, Table 3). This allows us to identify neutral Cu(II) complex as a form with the lowest sorption rate. Indeed, when we increased tartrate concentration to 0.01 M, and contribution of all Cu(II) anionic forms increased from 11.3 to 43.4% (Table 2), effective dynamic sorption capacity was twice as high as in 0.00312 M tartrate (3.87 over 1.73 mmol/g), even at a higher flow rate of 130 BV/h (Table 3). These data are not in agreement with conclusion made for Cu(II) sorption from citrate solutions on polyaminated resin about the higher sorption rate constant for the neutral CuHL 0 complex over anionic forms [16]. Regretfully, experimental data in [16] were obtained in batch and were not verified in fixed-bed models at different flow rates and metal/ligand ratios that could prove one or another hypothesis.
At the same time, and similarly to the works [16,18,19,46] on Cu(II) sorption on polyamine resins in the presence of citrate, we observed a drastic increase in maximal sorption capacity for Cu(II) in the presence of a low excess of tartrate ( Table 3). The maximal value of 5.61 mmol/g is in very good agreement with the sorption capacity of 5.5 mmol/g toward anionic dye earlier determined for this cryogel [35]. It should be mentioned that homogeneous distribution of Cu(II) in monolith cryogel after the sorption from both acetate and tartrate solutions was confirmed by SEM ( Figure S1, Supplementary information). Thus, the difference in sorption capacities and rate-dependence of the shape of the breakthrough curves can be attributed to the different mechanisms of Cu(II) sorption in the presence of acetate and tartrate.

Comments on Cu(II) Sorption Mechanism in Relation to Application of RCD-Complex Model
There are numerous examples in the literature of when sorption mechanism changes depending on the ionic form of the adsorbate, even at the same pH value. Acetate buffers with high concentrations are often used as a background media to investigate metal sorption properties, since it is believed that acetate does not affect metal binding due to the low stability of acetate complexes [18]. Although we have demonstrated above that acetate can have profound effect on Cu(II) sorption kinetics, a comparison of FT-IR spectra (Figure 3

Comments on Cu(II) Sorption Mechanism in Relation to Application of RCD-Complex Model
There are numerous examples in the literature of when sorption mechanism changes depending on the ionic form of the adsorbate, even at the same pH value. Acetate buffers with high concentrations are often used as a background media to investigate metal sorption properties, since it is believed that acetate does not affect metal binding due to the low stability of acetate complexes [18]. Although we have demonstrated above that acetate can have profound effect on Cu(II) sorption kinetics, a comparison of FT-IR spectra (Figure 3) of PEI cryogels after Cu(II) sorption from water and acetate solutions reveals the following similar features: band shifts of N-H bond bending (1300-1315 and 1560 cm −1 ) and C-N bond stretching (1100-1050 cm −1 ) Cu(II) chelation by a different type of PEI amino groups. The notable difference was observed only in the region around 1650 cm −1 , where bands of N-H bending can overlap with stretching vibration of carboxylic group. The widening and asymmetry of this band in PEI FT-IR spectra after Cu(II) sorption from both 0.1 M and 1 M acetate solutions does not exclude possibility of replacement of one water molecule with acetate in the coordination sphere of the PEI-Cu(II) complex, where copper has coordination number 5 [47]. However, it is important to emphasize that FT-IR spectra of PEI after the Cu(II) sorption in presence of acetate have the same features regardless of sorption conditions (ligand concentration and flow rate). In contrast, FT-IR spectra of PEI after Cu(II) sorption from 0.1 M, 0.01 M, and 3.12 mM tartrate solutions (spectra 6-8, Figure 3) feature significant changes in relative intensities and positions of the bands at 1315-1345 cm −1 and 1580-1650 cm −1 regions that can be interpreted in terms of overlapping N-H bending with asymmetric and symmetric stretching of carboxyl groups. Moreover, the shape of the peaks and area ratios depend on sorption conditions ( Figure S2, Supplementary materials) that can reflect changes in the binding mode depending on ionic speciation [48]. Cu(II) sorption from tartrate solutions can be additionally complicated by formation of polynuclear complexes. Calculations of Cu(II) ionic speciation in 3 mM tartrate solution with a tartrate:Cu(II) molar ratio of 2:1: using stability constants for bi-nuclear tartrate complexes [49] showed that 7.45% of Cu(II) In contrast, FT-IR spectra of PEI after Cu(II) sorption from 0.1 M, 0.01 M, and 3.12 mM tartrate solutions (spectra 6-8, Figure 3) feature significant changes in relative intensities and positions of the bands at 1315-1345 cm −1 and 1580-1650 cm −1 regions that can be interpreted in terms of overlapping N-H bending with asymmetric and symmetric stretching of carboxyl groups. Moreover, the shape of the peaks and area ratios depend on sorption conditions ( Figure S2, Supplementary materials) that can reflect changes in the binding mode depending on ionic speciation [48]. Cu(II) sorption from tartrate solutions can be additionally complicated by formation of polynuclear complexes. Calculations of Cu(II) ionic speciation in 3 mM tartrate solution with a tartrate:Cu(II) molar ratio of 2:1: using stability constants for bi-nuclear tartrate complexes [49] showed that 7.45% of Cu(II) exists in the [Cu 2 L 2 ] 0 form (Table S3, Supplementary materials). Formation of such complexes can contribute to the increase in sorption capacity in solutions with low tartrate concentrations.
However, Ling et al. [16] suggested that neutral Cu(II)-citrate complex was adsorbed on polyaminated resin with the release of the ligand via coordination with amine sites. The FT-IR spectrum of PEI after adsorption of Cu(II) in the presence of 3.12 mM and 10 mM of tartrates, containing 46.5% and 41.9% of cationic and neutral Cu forms, significantly differs from the PEI-Cu(II) spectrum. Thus, despite relatively low stability constants of Cu(II)-tartrate complex in comparison with PEI-Cu(II) complex, the mechanism of Cu(II) sorption from tartrate solution involves the formation of mixed ligand complexes, whose composition depends on sorption conditions. This fact can limit applicability of the RCDcomplex model to a Cu(II) tartrate system, since we hypothesized that one can modify RCD function obtained in water for predictive modeling of sorption in the presence of ligands if the type of Cu(II) binding to the sorption site remains the same.

Description of PEI Sorption Site Characteristics in Presence of Complexing Agents Using RCD Model
First, maps of PEI binding cites distribution in the space of constants of sorption (Ks) and desorption (Kd) rates were calculated using RCD model from breakthrough curves of Cu(II) sorption from water, 0.1 M acetate and 0.1 M tartrate. Under these sorption conditions, RCD function will be calculated for the sorption of [CuL 4 ] 6− form (>98%) in tartrate solution and mixture of four Cu(II) ionic forms (Table 1) in acetate solution. Site distribution patterns (Figure 4a) with at least two main site types, which differ in sorption rate and affinity, are similar for sorption from all three solutions. This indicates that the same functional fragments were involved in Cu(II) binding. However, one observes a gradual shift of the pattern to the lower sorption rate constants in the water-tartrate-acetate row. Only a small population of the "fast" sorption sites (Ks values > −1) was identified in the acetate solution. Despite an obvious difference in sorption kinetics, Q max values for Cu(II) sorption in 0.1 M acetate and water calculated using Equation (6) of the RCD model ( Figure 2b) do not differ significantly and corroborate with maximal dynamic sorption capacities (Table 3). Theoretical isotherms in 0.1 M acetate and tartrate solutions are also in good agreement with experimental data obtained in batch for equilibration time of 72 h (Figure 4b). This shows that all Cu(II) ionic forms in can be adsorbed by PEI but with different rates, and the rate-dependence is more pronounced in 0.1 M acetate solution.
Participation of the same functional groups in sorption of Cu(II) from solutions with and without complexing agents and the same range of affinity for the sorption centers under these sorption conditions (Figure 4a) allowed for the assumption that RCD function for Cu(II) sorption from water can be used to describe sorption in the presence of complexing agents and predict shape of the breakthrough curves at different flow rates and metal/ligand concentrations, if coefficients a i,j . in Equation (13) are found. However, a difference in sorption mechanism of Cu(II) from water and tartrate can negatively affect possibility to calculate one RCD function for the broad range of ligand concentrations, as one can see that the dynamic sorption capacities for Cu(II) significantly depend on the excess of tartrate (Table 3).

Determination of Sorption Rate Constants for Different Ionic Forms of Cu(II)
Although site distribution maps for Cu(II) sorption in the absence and presence of complexing agents (Figure 4a) do not differ drastically, the direct application of RCD function for Cu(II) adsorption from water for predictive modeling of Cu(II) sorption breakthrough curves in acetate and tartrate solutions does not give satisfactory results ( Figure 5). This is especially the case for Cu(II) sorption from acetate solution (Figure 5a) when RCD function for sorption in water cannot predict even the general trend of a breakthrough curve shape and position depending on the flow rate.

Determination of Sorption Rate Constants for Different Ionic Forms of Cu(II)
Although site distribution maps for Cu(II) sorption in the absence and presence of complexing agents (Figure 4a) do not differ drastically, the direct application of RCD function for Cu(II) adsorption from water for predictive modeling of Cu(II) sorption breakthrough curves in acetate and tartrate solutions does not give satisfactory results ( Figure  5). This is especially the case for Cu(II) sorption from acetate solution (Figure 5a) when RCD function for sorption in water cannot predict even the general trend of a breakthrough curve shape and position depending on the flow rate.

Determination of Sorption Rate Constants for Different Ionic Forms of Cu(II)
Although site distribution maps for Cu(II) sorption in the absence and presence of complexing agents (Figure 4a) do not differ drastically, the direct application of RCD function for Cu(II) adsorption from water for predictive modeling of Cu(II) sorption breakthrough curves in acetate and tartrate solutions does not give satisfactory results ( Figure  5). This is especially the case for Cu(II) sorption from acetate solution (Figure 5a) when RCD function for sorption in water cannot predict even the general trend of a breakthrough curve shape and position depending on the flow rate.   (Table 1), shows a strong dependence on the flow rate. This suggests significant differences in sorption rate constants for these various ionic forms. Additional experimental data were obtained for the sorption at Cu(II) excess (0.001 M acetate) and high acetate excess (1 M acetate) to process more breakthrough curves simultaneously using RCD model and obtain one set of a ij coefficients in a Tylor series (Equation (13)), which would fit all breakthrough curves obtained in solutions with different ratios of Cu ionic forms. The Tylor series coefficients calculated for this data set for modifying RCD function for Cu(II) sorption from water are given in Table S1 (Supplementary Materials). It should be mentioned that to shorten computation time, Kd values were limited to −2.5. This is reasonable assumption, since Figure 3a shows that the majority of PEI sorption sites under the chosen sorption conditions fall within this range of desorption rate constants. Figure 6a shows that the RCD-complex model provides a good fit for all breakthrough curves of Cu(II) sorption from acetate solutions varying in ligand concentrations. The corresponding distribution of different Cu(II) ionic forms in the space of sorption/desorption rate constants (Figure 6b) shows that sorption rate constants of the Cu 2+ ions are notably higher in comparison with those of other ionic forms. The most slowly adsorbed form is a neutral complex [CuL 2 ] 0 that can be attributed to its lack of electrostatic attraction to the amino groups of PEI, which are positively charged at pH 5.2. ferences in sorption rate constants for these various ionic forms. Additional experimental data were obtained for the sorption at Cu(II) excess (0.001 M acetate) and high acetate excess (1 M acetate) to process more breakthrough curves simultaneously using RCD model and obtain one set of aij coefficients in a Tylor series (Equation (13)), which would fit all breakthrough curves obtained in solutions with different ratios of Cu ionic forms. The Tylor series coefficients calculated for this data set for modifying RCD function for Cu(II) sorption from water are given in Table S1 (Supplementary Materials). It should be mentioned that to shorten computation time, Kd values were limited to −2.5. This is reasonable assumption, since Figure 3a shows that the majority of PEI sorption sites under the chosen sorption conditions fall within this range of desorption rate constants. Figure 6a shows that the RCD-complex model provides a good fit for all breakthrough curves of Cu(II) sorption from acetate solutions varying in ligand concentrations. The corresponding distribution of different Cu(II) ionic forms in the space of sorption/desorption rate constants (Figure 6b) shows that sorption rate constants of the Cu 2+ ions are notably higher in comparison with those of other ionic forms. The most slowly adsorbed form is a neutral complex [CuL2] 0 that can be attributed to its lack of electrostatic attraction to the amino groups of PEI, which are positively charged at pH 5.2. To prove the prognostic value of the RCD-complex model, we have simulated a breakthrough curve of Cu(II) sorption from the solution with Cu(II) and acetate concentrations, which were not used in the experiment designed for the calculation of RCDcomplex function. The result of predictive breakthrough curve modeling using the RCD function for Cu(II) sorption from water and coefficients of Tylor row is shown in Figure 7a, which demonstrates very good agreement between experimental and theoretical data.
When we extended range of the ligand concentration in tartrate solutions, we failed to find one set of the Tylor row coefficients, which would fit all experimental curves (Figure 8a). At a low tartrate concentration, corresponding to a 2:1 tartrate: Cu(II) mole ratio (Cu speciation is given in Table 2), the breakthrough curve had a very low slope after the To prove the prognostic value of the RCD-complex model, we have simulated a breakthrough curve of Cu(II) sorption from the solution with Cu(II) and acetate concentrations, which were not used in the experiment designed for the calculation of RCD complex function. The result of predictive breakthrough curve modeling using the RCD function for Cu(II) sorption from water and coefficients of Tylor row is shown in Figure 7a, which demonstrates very good agreement between experimental and theoretical data.
When we extended range of the ligand concentration in tartrate solutions, we failed to find one set of the Tylor row coefficients, which would fit all experimental curves (Figure 8a). At a low tartrate concentration, corresponding to a 2:1 tartrate: Cu(II) mole ratio (Cu speciation is given in Table 2), the breakthrough curve had a very low slope after the breakpoint, while dynamic sorption capacity increased drastically in comparison with sorption from water and 0.1 M tartrate (Table 3). Thus, we have calculated the Tylor row coefficients (Table S2, Supplementary Materials) and distribution of different ion forms in the space of sorption/desorption rate constants for the set of experimental curves obtained in 0.1 M tartrate only (Figure 8b). breakpoint, while dynamic sorption capacity increased drastically in comparison with sorption from water and 0.1 M tartrate (Table 3). Thus, we have calculated the Tylor row coefficients (Table S2, Supplementary Materials) and distribution of different ion forms in the space of sorption/desorption rate constants for the set of experimental curves obtained in 0.1 M tartrate only (Figure 8b). It can be observed that the difference between the sorption rate constant of the free metal ion and other ionic forms, except for [CuL2] 2-, is less pronounced than in Cu(II)acetate solutions. This observation may indicate a disparity in the Cu(II) sorption mechanism in the presence of these two ligands, which has been demonstrated with FT-IR spectroscopy ( Figure 3, Section 2.4). Expectedly, the RCD-complex model with the set of parameters obtained in 0.1 M tartrate failed to predict the shape of the breakthrough curve at low tartrate concentration (Figure 7b).  It can be observed that the difference between the sorption rate constant of the free metal ion and other ionic forms, except for [CuL 2 ] 2− , is less pronounced than in Cu(II)acetate solutions. This observation may indicate a disparity in the Cu(II) sorption mechanism in the presence of these two ligands, which has been demonstrated with FT-IR spectroscopy (Figure 3, Section 2.4). Expectedly, the RCD-complex model with the set of parameters obtained in 0.1 M tartrate failed to predict the shape of the breakthrough curve at low tartrate concentration (Figure 7b).
At this point, it is important to emphasize that assumptions made in the RCD-complex model, first of all, are based on the fact that despite the difference in ionic speciation of the metal ion, the adsorbate binds to the sorbent via the same mechanisms in water and in the presence of a complexing agent. In other words, sorption of metal ion and metal-chelate results in binding of a metal ion by the same functional groups of the sorbent and the release of the ligand in the latter case. Most likely, this necessary condition is not fulfilled in a broad range of tartrate concentrations. In such cases, the necessary conditions for application of the RCD-complex model are not fulfilled, making it impossible to find a single set of parameters for predictive modeling of sorption kinetics.
Thus, the decision on applicability of RCD-complex model for predictive modelling of sorption kinetics in a broad range of solution compositions and sorption conditions should be based on information about the sorption mechanism. If all ionic forms of metal ion form with the sorption center complexes of the same or a very similar configuration, the RCD-complex model can be used for predictive modeling of sorption dynamics.

Cryogel Synthesis and Characterization
The PEI cryogels were obtained as described in detail in [35]. Briefly, the solution for cryogel fabrication was prepared by addition of the cross-linker at the DGE-1,4-BD:PEI molar ratio 1:4 to the 5% PEI solution under intensive stirring at room temperature. The plastic syringes of diameters (0.47 cm) were filled with the mixed solution and kept frozen at −20 • C for 7 days. After thawing at room temperature, the obtained PEI cryogels were washed with 10 mL of 0.1 M NaOH solution and distilled water using peristaltic pump (Ismatec, Wertheim, Germany).
Fourier transform infrared (FT-IR) spectra for air-dried PEI cryogel before and after Cu(II) sorption were recorded using an IR Affinity-1 spectrometer with QATR 10 single reflection ATR accessory (Shimadzu, Kyoto, Japan). The monolith cryogels after the Cu(II) sorption were thoroughly washed with a large volume of distilled water to remove all unbound ions.

Investigation of the Sorption in Fixed-Bed and Sorption Isotherms
The dynamics of Cu(II) ion sorption on monolith PEI cryogel was investigated by pumping solution of copper nitrate in water, sodium acetate, or potassium-sodium tartrate solutions through the sorption column at different flow rates. Prior to solution feeding, the column was equilibrated by passing 10 mL of the background solution without Cu(II) added. The samples were collected every 5 mL, and the copper concentrations were determined by atomic absorption spectrometry (AAS) using a AA-6200 Atomic Absorption Flame Emission Spectrophotometer (Shimadzu, Kyoto, Japan) device.
The sorption isotherms were investigated at solid:liquid ratio 1:1000, the contact time was 72 h, and the solutions with sorbents were agitated using a Biosan PSU-20i orbital shaker (Riga, Latvia) at 200 rpm. At least three replicates were made to ensure the results reproducibility. The adsorbed amounts were calculated using the difference in initial and equilibrium concentrations of the metal ions in the solutions determined by AAS.

Calculation of RCD Functions for Metal Ions in the Presence of Complexing Agents
To obtain the RCD function from sorption dynamics experiments, several breakthrough curves obtained at different flow rates and fixed ligand concentration (0.1 M) were processed using the same RCD model, which was used to describe Cu(II) sorption in a fixed-bed application in the absence of complexing agent [29]. Using simulation modeling as described in [29], an averaged 3D distribution functions for Cu(II) sorption sites of PEI in presence of complexing agents (acetate and tartrate) were plotted as contour lines in the space of sorption/desorption rates constants-ρ(pK s, p K d ),pK s , pK d .
Alternatively, experimental breakthrough curves obtained at different ligands concentrations and flow rates were processed using RCD function obtained earlier for Cu(II) sorption from aqueous solution [29], and Equation (15) to calculate the coefficients a ij for every complex form. On figures describing the distribution of sorption centers of different types in the pK s -pK d space, the deter position coincides with the circle center. Here, the circle radius is proportional to the maximal capacity of this very center (q max ).

Predictive Modeling Using RCD Model
The predictive modeling of the breakthrough curves for different initial concentrations of the adsorbate and ligands have been performed using the RCD function determined from fixed-bed experiments for Cu(II) from water or using RCD complex function calculated from fixed-bed experimental data for solutions at different ligand concentrations and flow rates.

Conclusions
Aside from affinity, selectivity, and high-sorption capacity, the kinetic characteristics of sorbents are of utmost importance, as they determine productivity of both industrial water treatment facilities and small-size, point-of-use (POU) filters operated at high flow rates. In this work, we have discussed how and when one can predict breakthrough curves of metal ion sorption in the presence of complexing agents, if diffusion limitations are eliminated via well-developed porous structure of the sorbent, and if chemisorption is a rate-limiting stage.
A new version of the earlier developed and verified extended Rate Constants Distribution (RCD) model [27][28][29], which additionally includes complex-formation equilibria, was discussed here and applied to an investigation of Cu(II) sorption kinetics on a supermacroporous monolith polyethyleneimine (PEI) cryogel in the presence of acetate and tartrate. The RCD-complex model assumes that in order to predict the shape of the breakthrough curves in the presence of a ligand, one can use the parameters (RCD-function) determined for the metal ion sorption from ligand-free solutions along with a set of coefficients, which account for the variation in sorption rate constants of different metal ionic forms. However, this approach assumes the same or a very similar configuration of the metal-sorbent site after adsorption from solution with and without a complexing agent. Thus, information on sorption mechanism is crucial for the decision on applicability of the model.
The application of the RCD-complex model to the sets of Cu(II) sorption breakthrough curves obtained at different flow rates and metal speciation in tartrate and acetate solutions allowed us to conclude the following points: (i) despite the lower stability of Cu(II)-acetate complex, it can have a more pronounced negative effect on sorption kinetics than tartrate at high ligand:metal ratios; (ii) the neutral Cu(II)-acetate complex adsorbs with the lowest rate; (iii) the presence of tartrate results in a twofold increase in the dynamic sorption capacity of PEI for Cu(II) at low ligand:metal ratios even at high flow rates (up to 130 BV/h); and (iv) the RCD-complex model successfully predicts breakthrough curves of Cu(II) sorption in the broad concentration range of acetate but not of tartrate. The latter fact was related to alteration of sorption mechanism in the presence of tartrate and formation of the mixed ligand Cu(II)-PEI-tartrate complexes instead of Cu(II)-PEI complex, despite a large difference between formation constants of Cu(II)-PEI and Cu(II)-tartrate complexes.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Expansion in a Taylor series is conducted around the point pKs 0 Me , pKd 0 Me : After removal of brackets and grouping of terms by equal powers in (A3), one obtains: pKs MeL n = a n 00 + a n 10 pKs Me + a n 01 pKd Me + a n 20 pKs Me 2 + a n 11 pKs Me pKd Me + . . .